Optical conductivity in cluster dynamical mean field theory: formalism and 
application to high temperature superconductors 
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The optical conductivity of the one-band Hubbard model is calculated using the 'Dynamical 
Cluster Approximation' implementation of dynamical mean field theory for parameters appropriate 
to high temperature copper-oxide superconductors. The calculation includes vertex corrections and 
the result demonstrates their importance. At densities of one electron per site, an insulating state 
is found with gap value and above-gap absorption consistent with measurements. As carriers are 
added the above gap conductivity rapidly weakens and a three component structure emerges, with 
a low frequency 'Drude' peak, a mid-infrared absorption, and a remnant of the insulating gap. The 
mid-infrared feature obtained at intermediate dopings is shown to arise from a pseudogap structure 
in the density of states. On further doping the conductivity evolves to the Drude peak plus weakly 
frequency dependent tail structure expected for less strongly correlated metals. 

PACS numbers; 74.25.Gz, 72.80.-r, 74.25. Fy, 71.27.-l-a 
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The frequency-dependent ('opticar) conductivity a{Q) 
is an important probe of electronic condensed matter 
physics, revealing electronic band gaps, scattering pro- 
cesses, carrier number and effective mass. Optical con- 
ductivity measurements have played a particularly im- 
portant role in the study of high temperature copper- 
oxide superconductorsii^, revealing behavior which de- 
viates sharply from conventional band theoretic expec- 
tations. In the undoped 'parent compounds' such as 
La2CuOi measurements^ reveal a large ~ 1.75eV gap 
which persists essentially unchanged as temperature is 
raised above the Neel temperature ~ 300i^ while 
band theory predicts metallic behavior in the absence 
of antiferromagnetism. As the materials arc doped, the 
gap feature weakens and absorption appears at lower fre- 
quencies. The lower frequency absorption is often^ de- 
composed into two parts, a 'Drude' peak centered at 
i7 = and a 'mid-IR' structure at fl 0.5eV; both 
parts are characterized by an oscillator strength (integral 
of conductivity over a frequency range) which is small 
compared to the band theory value and increases as the 
doping is increased^. The 'mid-IR' band has been var- 
iously interpreted as an effect of scattering of carriers 
from spin fluctuations^!^ or other bosons^i^, a signature 
of two-component absorption^ and as an indication of a 
novel charge 2e excitatio n^°d^ , but a clear consensus on 
the interpretation has not emerged. 

Cluster dynamical mean field approximation a^^-^^ to 
the one-band Hubbard model have been argued to pro- 
vide a reasonable description of the physics of the 
cuprates. Single electron properties such as the pho- 
toemission spectra have been argue d^^i^^'^^i^^'^^ to be 
in good agreement with data. However, in cluster dy- 
namical mean field theory evaluation of two particle re- 
sponse functions such the optical conductivity requires 
computation of a vertex function. In this paper we 
show that the vertex correction may be computed and 
makes a significant contribution (especially to the con- 
ductivity of undoped and lightly doped materials). Our 



— x=0 

— x=0.008 

— - x=0.054 
■ - x=0.089 

— x=0.128 

— x=0.21 




FIG. 1: Optical conductivity calculated for indicated dopings 
using 4-site DCA approximation to the paramagnetic phase of 
the two-dimensional one-band Hubbard model at parameters 
corresponding to high temperature superconductors at tem- 
perature T ~ 400K and other parameters described in the 
text. Inset: 2/tt times integral of optical conductivity over 
low frequency fl < 0.6t ('Drude', solid line, black on-line), in- 
termediate frequency f2 = 2t ('Mid-IR', long-dashed line, red 
on-line), and all frequencies ('Total', light solid line, green 
on-line), along with independently computed kinetic energy 
('KE', light dashed line, blue on-line). 



main result, a computation of the variation with dop- 
ing of the frequency dependent conductivity cr(f2) of 
the Hubbard model at physically relevant parameters, 
is shown in Fig. [TJ The curves bear a striking similar- 
ity to the conductivity of hole-doped high temperature 
superconductorsii^. The conductivity calculated for the 
undoped system displays a characteristic insulating spec- 
trum with a gap ~ 5t « l.SeV similar in magnitude 
to that observed in La2CuOi and an above-gap absorp- 



tion strength corresponding to ~ ISOOfJ 



about 



30% larger than observed. On doping, this gap is rapidly 
destabilized, again in a manner consistent with measure- 
ments. Absorption in the near-gap region is suppressed, 
while substantial absorption strength appears at low fre- 
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quencies. As the doping is increased the optical response 
may be described in terms of a 'Drude' peak indicat- 
ing coherent quasiparticlc motion, an additional 'mid- 
IR' feature {lo ~ 0.5t) and a high frequency tail. By 
the highest doping the conductivity has evolved to the 
'Drude peak' plus weak high frequency tail characteristic 
of Fermi-liquid metals. 

To further characterize the evolution of the conduc- 
tivity we present in the inset of Fig. [1] the partial opti- 
cal integrals K{Vl) = ^ duja{uj) for a low frequency 
(fl = 0.6t encompassing the Drude peak), an intermedi- 
ate frequency ^ 2t , somewhat less than half of the 
gap value and the total integral 51 — > oo (which is seen to 
agree with the independently calculated kinetic energy). 
The shift of spectral weight, first into a mid-IR band and 
then into a Drude peak is very similar to observations in 
high- Tc cuprates, although the calculated mid-IR spec- 
tral weight is higher than values inferred from data (see 
e.g. Fig. 1 in Refs. Uli.) 

The derivation of our results begins from the current- 
current response function Xjjit) = idt(j{t) which relates a 
spatially uniform, time dependent (therefore transverse) 
electric field E to the current j it creates. For a system 
described by a Hamiltonian H = T + U with interac- 
tions U which depend only on particle and spin densi- 
ties (not on particle or spin currents) we have j{t) = 
Tr[J(A)G(t; {A})]. The current operator J is obtained 
from the derivative of the single-particle terms T with re- 
spect to vector potential: J = ST /5A while the electron 
Green function G(t; {A]) = {idt - T({A}) - S({l}))-i 
is to be computed in the presence of the time depen- 
dent vector potential A. Here bold face quantities de- 
note matrices in the space of states of the system (in- 
cluding the spatial indices). Expanding to linear or- 
der in A and introducing the 'kinetic energy' operator 
K = 5'^T/5A^ and vertex operator T = 6Y,/5A we ob- 
tain Xj] = Xdia + Xbubblc + Xvertox with 

Xdia(i - t') = Tr [KG(t = 0)] 5{t - t'), (1) 
Xbubbie(i - t') = Tr [jG(t - t')3G{t' - t)] , (2) 

Xvertex(i " t') = Tr [jG(t - h)f{ti - t\ t' - t2)G{t2 - t) 

(3) ■ 

(convolution on internal time indices is to be under- 
stood). Combining the Kramcrs-Kronig relation between 
real ( ) and imaginary ( ) parts of cr with the gaugc- 
invariance condition that for a non-superconducting ma- 
terial XjA^ = 0) = impliesi^ Xdia = / ^Xi'j,{<^)/uJ = 
/ ^cr [lo), which is the usual f-sum rule. 

In the dynamical cluster approximation (DCA)i^i^ 
implementation of the dynamical mean field 
approximation^! one tiles the Brillouin zone into 
a = 1...N non-overlapping equal- area regions and 
approximates the electron self-energy E(fc,a;) by the 



piecewise constant form 



N 



(4) 



with (j)a{k) = 1 if fc is in tile a and zero otherwise. Thus if 
N ^ I the self-energy has an explicit momentum depen- 
dence arising from the discontinuities at the boundaries 
of the tiles. The self energies Sa are computed from the 
solution of an A'^-site quantum impurity model which in- 
volves the interactions of the original model (projected 
onto the impurity cluster) and mean field functions Q^^ 
which are fixed by the self consistency equations 



5a ^ = Sa 



{dk)G{k) 



(5) 



where the integral is over the tile a with appropriate mea- 
sure [dk), G{k) is the Green function of the lattice prob- 
lem computed with S defined by Eq. 

For the conductivity we require the vertex function 
r(a; + Q,,uj) = 5Y./5A, which has two sources: the ex- 
plicit dependence on k arising from the momentum-space 
discontinuities and any additional dependence of S on 
A. The additional A-dependence arises via the impurity 
model from a dependence of Q^^ on A which may be 
computed by linearizing Eq. ([5]) in A. Denoting the first 
order changes in Q and S by 5Q and 5Ti we have 



5g-^~H{5iio.}]^~i':-A 

with (time arguments are not written explicitly) 



(6) 



t = ^j{dk)G{k)^G{k)^ 1 (7) 

/s = 5^'^ + G-'(^jmG{k)5^^G{k)^G-' (8) 

The canonical tiling for clusters of size A^ = 1,2,4 
produces momentum sectors with symmetry such that 
h=Q. Thus as noted by Ref. US (for A^ = 1) andlH (for 
A^ = 4), in these clusters there is no explicit dependence 
of Sa on A. Ref. [H further argued that for TV = 4 all 
vertex corrections vanished. This is incorrect, although 
the vertex corrections turn out to be unimportant for the 
situation of interest to Ref. |23| . The explicit momentum 
dependence provides a non-vanishing vertex correction 
arising from the momentum space discontinuities which 
occur along the lines k°'^ separating tile b (on the of larger 
k side) from tile a (on the smaller k side) . We define n"'' 
to be the normal to this line. In the = limit the ver- 
tex correction is directly given by dTi/dk. To determine 
the vertex correction for we consider the Ward 
identity QVp- q- fj = G-^{k + q,uj + n) - G-^{k,uj) 
relating the charge and current Fj vertices to the in- 
verse Green functions. Within the DCA approximation 
Tp arises from the functional derivative of S with respect 
to a time dependent chemical potential; it is computed 
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FIG. 2: Main panel: conductivity calculated from 4-site 
DCA approximation to one-band Hubbard model for U = 9t, 
i! = — 0.3t and dopings a; = (heavy lines, red on-line) and 
0.054 (light lines, black on-line) by continuing self energy. Full 
lines; conductivity; dashed lines; contribution to conductiv- 
ity from convolution of bubble diagrams Eq. [2] dash dotted 
lines; contribution from vertex corrections, Eq. O Larger In- 
set; comparison of conductivity computed by continuing self 
energy (dashed line, red on-line) and continuing Matsubara 
response function (solid line, black on-line) at f/ = 6t and 
t' — 0. Smaller inset; expansion of high frequency negative 
conductivity region for U = 6t and t' = 0. 

along the lines of Eq. © but with dsk/dk- A replaced by 
dfi. Because this perturbation is a scalar it has no contri- 
bution from the functions (f>k and therefore cannot have 
any term proportional to a delta function of k. Thus we 
may identify the current vertex from the contributions 
in the Ward identity proportional to delta functions in 
fc-space, yielding 

f'iLO+n, u) = n''^ {T.b{LU + n)- T.a{uj)) S ((fc - r'') ■ n"'') 

(9) 

We apply the formalism to the 4 site DCA ap- 
proximation to the two dimensional Hubbard model, 

—2t (cos fcj; + COS ky)~At' cos cos ky. The current oper- 
ator for the X direction is jx = 2t sin kx + 4t' sin kx cos ky 
and the kinetic energy operator K = 2tco&kx + 
4t' cos fca; cos /cy . For this approximation, Eq. [9] is the 
appropriate vertex correction. We restrict attention to 
the paramagnetic phase. We use the numerically ex- 
act continuous- time auxiliary fieldSi (CT-AUX) impu- 
rity solver to solve the impurity model and construct 
vertex functions and the conductivity. Parameters rel- 
evant to high temperature copper-oxide superconductors 
are t « 0.35el/, t' « -0.3*^5 and U - 9ii. The precision 
of calculations for t' ^ or n ^ 1 is limited by a fermion 
sign problem; we therefore present results for U =^ Qt and 
t' = Q where higher precision data can be obtained. 

Because we solve the model on the imaginary axis, an 
analytical continuation is required to obtained real fre- 
quency information. This may be done in two ways: ei- 
ther by continuing the self energies^ and then computing 
the conductivity or by continuing directly the Matsubara- 



axis conductivity. Analytical continuation requires a co- 
variance matrix of error estimates. For the self energy the 
covariance matrix is available from the QMC data, while 
for the conductivity we estimate the covariance (which is 
highly non-diagonal because of the 1/57 in the definition) 
from an ensemble of 16 independent QMC solutions of 
the mean field equations. The continuations arc tested by 
comparing the directly computed Matsubara axis cr(i57„) 
to the same quantity 'back-computed' from the contin- 
uation. Significant differences in real-axis conductivity 
correspond to variations of a few times 10~^ in (T(if2„), 
setting a stringent requirement on the quality of the data. 
The inset of Fig. [2] presents highly precise data obtained 
at the sign-problem- free parameters U = 6t and t' = 0; 
we see that the two continuation methods yield similar 
results; differences between them are a measure of the 
best-case uncertainties in the continued a. At U = 9t 
and t' = — 0.3t the differences between the methods are 
larger; in particular the gap edges are much broader in 
the traces obtained by continuing the Matsubara axis re- 
sponse functions; we believe this broadening is unphys- 
ical for the reasons given in Ref. H^. The conductivity 
obtained from continuing cr(if2„) is found in general to 
produce a back-continued a in worse agreement with orig- 
inal data than the conductivity obtained by continuing 
the self energy; we therefore present results obtained from 
the latter method. One difficulty must be noted. As can 
be seen from Fig. [2l in the high frequency regime (fre- 
quencies well above the gap) the vertex correction acts 
to cut off the high frequency tail found in Xbubbic, adding 
a negative contribution to the positive-definite Xbubbie so 
that the total contribution nearly vanishes. In calcula- 
tions based on continuing the self energy first, the vertex 
correction in fact overcompensates, leading to an unphys- 
ical negative conductivity for some high frequencies (if 
a{ifln) is directly continued the a (uj) is by construction 
positive) . We believe that the overcompensation is a nu- 
merical artifact. Our numerical uncertainties, both from 
the QMC measurement and the analytical continuation, 
are largest in this regime. The magnitude of the unphys- 
ical negative contribution is small: for our U = 9t calcu- 
lations the spectral weight in the negative region ranges 
from 5% of the total spectral weight at a; = to 0.4% 
at our highest doping. Further, the overshoot is found to 
decrease as the numerical accuracy of our computation is 
improved and as seen from the second inset to Fig. [5] can 
be made smaller than 2% of the total spectral weight in 
our best case. However, to date we have been unable to 
eliminate the negative region entirely. 

The main panel of Fig. [5] presents the conductivity for 
t' ~ — 0.3t and U = 9i as well as its decomposition into 
'bubble' (Eq. [2]) and 'vertex' (Eq. [3]) contributions for 
dopings a; = and x = 0.054. The vertex correction is 
seen to make a non-negligible contribution to the conduc- 
tivity and to be essential to fulfilling the /-sum rule. It 
decreases in importance as doping increases, as expected 
because with increased doping the self energy becomes 
more isotropic in momentum space. This aspect of our 
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FIG. 3: Main panel: Density of states calculated for the sector 
containing the Fermi surface from 4-site DCA approximation 
at U — 9t and t' = —0.3t at dopings indicated. Inset; expan- 
sion of the near-Fermi-surface region. 

result disagrees with Ref. [ll| which stated (on the basis 
of a comparison of a 4-site CDMFT conductivity com- 
puted without vertex corrections to the f-sum rule) that 
vertex corrections were unimportant near half filling and 
increased in importance as the doping increased; on the 
other hand at lower frequencies and higher dopings where 
vertex corrections are of less importance our result is sim- 
ilar to that of Ref. [ll|. The origin of the difference is not 
clear. 

At a; = the vertex corrections bring two effects: they 
increase the magnitude of the conductivity in the above- 
gap region and they steepen the rate at which the con- 
ductivity rises above the gap edge. We believe that these 
two effects are consequences of the short ranged order 
captured by the DCA approximation. Inspection of the 
Green function (not shown here) indicates that on the 
single-particle level the state is an indirect-gap insulator 
with the highest energy filled states at a different momen- 
tum from the lowest energy empty states. As shown e.g. 
in Ref. and [H, the backfolding associated with long- 
ranged order converts the indirect gap to a direct one, 
dramatically steepening the conductivity onset. The ver- 
tex correction provides a similar effect. The vertex cor- 
rection also expresses the 'coherence factor' physics as- 
sociated with strong two-sublattice spatial correlations. 
To see this, consider a mean-field model of a material 
with two-sublattice order, for which the electron propa- 
gator has 'normal' (G ^ (cfccj^)) and 'anomalous' [F ~ 



(cfc4-KQ)) parts given by G{k,uj) = [uj ~ ek+Q)/{{uj - 
ek)(uj - Ek+q) - A^) = \/{uj -Sk- A2/(a; - ek+q)) and 
F{k,Lu) = AV((w-efe)(w-efc+Q)-A2) respectively. The 
conductivity in the ordered state is computed from the 
sum of a 'G — G' and an 'F — bubble which give equal 
contributions to the conductivity for frequencies near the 
gap edge. The state uncovered in the 4-site DCA calcula- 
tion has no long ranged order, so the anomalous [F) part 
vanishes and convolution of bubble diagrams would cap- 
ture only the G ~G contribution. The vertex corrections 
in effect add back the F ~ F term. 

The conductivity is related to the electron spectral 
function, shown in Fig. O The initial doping moves 
the chemical potential into the lower Hubbard band and 
rapidly broadens the sharp peak at the edge of the up- 
per Hubbard band; we see also from Fig.[T]that the form 
of the near-gap-edge conductivity changes substantially. 
Interestingly, in the doped materials the remains of the 
above-gap absorption is entirely expressed by the vertex 
correction. At intermediate dopings x = 0.054 and 0.089 
the many-body density of states exhibits a 'pseudogap', a 
small gap at the Fermi level previously notedi^iiii^l. Ex- 
citations across the pseudogap have the correct energy to 
account for the mid-IR feature observed in the data and 
in high-Te materials (a similar connection was made in 

Rcf[iaini). 

To summarize, we have presented theoretically consis- 
tent calculations of the optical conductivity of the Hub- 
bard model within the 'DCA' implementation of cluster 
dynamical mean field theory. The calculated results bear 
a very great similarity to the conductivity observed in 
high temperature superconductors. The important role 
played by spatial correlations in the cluster DMFT ap- 
proximation (expressed in the calculation by vertex cor- 
rections) is seen from the rapid rise of the conductivity 
above the gap edge and the rapid changes with doping, 
while vertex corrections are less important at higher dop- 
ing and lower frequency. Important directions for future 
study include extensions to the case of Raman scattering 
and to larger clusters. 
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